<html>
<!--
    This work is licensed under Creative Commons GNU LGPL License.

    License: http://creativecommons.org/licenses/LGPL/2.1/
    Version: 0.5
    Author:  Michael Shook/2009
    Web:     http://del.icio.us/mshook
-->
<head>
<title>Filter Design - Calculate and display the impulse response (kernel) and frequency response for a windowed-sinc FIR filter - firkernel.htm</title>
<style type="text/css">
body {
  font-family: sans-serif;
  font-size: 16px;
  margin: 50px;
}
</style>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="description" content="Filter Design - Calculate and display the impulse response (kernel) and frequency response for a windowed-sinc FIR filter - firkernel.htm">
<meta name="keywords" content="filter,javascript,digital,high pass,lowpass,band pass,FIR,dsp,signal,processing,signal processing,finite impulse response,finite,impulse,response,impulse response,frequency,frequency response,
band reject">
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <!--[if IE]><script language="javascript" type="text/javascript" src="../excanvas.pack.js"></script><![endif]-->
    <script language="javascript" type="text/javascript" src="jquery.js"></script>
    <script language="javascript" type="text/javascript" src="jquery.flot.js"></script>

<script>
// To do:
//  * Turn it into a Google Gadget?
FIR = {
    // All the filter types start with the low-pass filter.
    // The Basic code that I started with is here:
    //   http://www.dspguide.com/ch16/4.htm
    // Basically this 
    //  1. does an Inverse Discrete Fourier Transform (IDFT)
    //      on the desired frequency response (which is just a step
    //      function from 1 to 0 at the cutoff frequency - the result
    //      has the form of the sinc function: sin(x)/x) and
    //  2. apply the Hamming window function to smooth
    //      out the response, since we cut off the infinite impulse
    //      response.
    // An explanation of the sinc function is here:
    //   http://www.dspguide.com/ch11/2.htm
    // An explanation of the IDFT is here:
    //   Synthesis, Calculating the Inverse DFT
    //   http://www.dspguide.com/ch8/5.htm
    calculateLowpass: function(sampleRate, filterOrder, cutoffFrequency) {
        var cutoffRadians = 2 * Math.PI * cutoffFrequency / sampleRate;
        var sum = 0;
        var h = [];    // h[] will be the impulse response or kernel of the filter
        for (i=0; i<=filterOrder; i++) {
            if ( i-filterOrder/2 == 0 ) { 
                h[i] = cutoffRadians;
            }
            else {
                // This is the basic sinc function calculation
                h[i] = Math.sin(cutoffRadians * (i-filterOrder/2)) / (i-filterOrder/2);
                // And here we apply the Hamming window
                // The explanation is here:
                // http://www.dspguide.com/ch16/1.htm
                h[i] *= (0.54 - 0.46 * Math.cos(2*Math.PI*i/filterOrder));
            };
        };
        // Normalize the filter kernel for unity gain at DC
        for (i=0; i<=filterOrder; i++) {
            sum = sum + h[i];
        };
        for (i=0; i<=filterOrder; i++) {
            h[i] /= sum;
        };
        return h;
    },
    // See the bottom of
    // http://www.dspguide.com/ch14/4.htm
    // for an explanation of spectral inversion
    invertSpectrum: function(h) {
        var i;
        for (i=0; i<h.length; i++) {
             h[i] = -h[i];
        }
        h[(h.length-1)/2] += 1.0;
        return h;
    },
    // See:
    //   Chapter 8: The Discrete Fourier Transform
    //   http://www.dspguide.com/ch8/6.htm
    // for a very thorough explanation of the DFT.
    // The chapter also explains the Inverse DFT
    // (IDFT):
    // The Basic code I started with for the DFT is
    // here:
    //   Analysis, Calculating the DFT
    //   http://www.dspguide.com/ch8/6.htm
    DFT: function(signal) {
        var rex = []; // real part of the frequency domain
        var imx = []; // imaginary part of the frequency domain
        var mag = [];
        var phase = [];
        var n = signal.length;
        var m = n/2;
        var i = 0;
        var k = 0;
        for (k=0; k<=m; k++) {
            rex[k] = 0;
            imx[k] = 0;
        };
        for (k=0; k<=m; k++) {
            for (i=0; i<n; i++) {
                rex[k] += signal[i] * Math.cos(2*Math.PI*k*i/n);
                imx[k] -= signal[i] * Math.sin(2*Math.PI*k*i/n);
            };
        };
        // Convert to polar notation - See:
        //   http://www.dspguide.com/ch8/8.htm
        for (k=0; k<=m; k++) {
            mag[k] = Math.sqrt(rex[k]*rex[k]+imx[k]*imx[k]);
            if (rex[k] == 0.0) rex[k] = 1e-20;
            phase[k] = Math.atan2(rex[k],imx[k]);
        };
        return {"rex": rex, "imx": imx, "mag": mag, "phase": phase};
    },
    calculateHighpass: function(sampleRate, filterOrder, cutoffFrequency) {
        return FIR.invertSpectrum(FIR.calculateLowpass(sampleRate,
            filterOrder, cutoffFrequency));
    },
    calculateBandreject: function(sampleRate, filterOrder, cutoff1, cutoff2) {
        return FIR.addArrays(
            FIR.calculateLowpass(sampleRate, filterOrder, cutoff1),
            FIR.calculateHighpass(sampleRate, filterOrder, cutoff2));
    },
    calculateBandpass: function(sampleRate, filterOrder, cutoff1, cutoff2) {
        return FIR.invertSpectrum(FIR.calculateBandreject(sampleRate, filterOrder,
            cutoff1, cutoff2));
    },
    addArrays: function(a1, a2) {
        var i;
        var a = [];
        for (i=0; i<a1.length; i++) {
            a[i] = a1[i] + a2[i];
        };
        return a;
    },
    arrayOfZeros: function(n) {
        var i;
        var a = [];
        for (i=0; i<n; i++) {
            a[i] = 0.0;
        };
        return a;

    },
    formatBasicDSP: function(samplerate, filterorder, 
            cutoff1, cutoff2, filtertype, h) {
        var resultString = "";
        var i;
        resultString += "samplerate=" + samplerate + "\n";
        resultString += "# samplerate: " + samplerate;        resultString += " " + FIR.filtertype;
        resultString += " order: " + filterorder;
        resultString += " cutoffs: " + cutoff1 + " " + cutoff2 +"\n";
        resultString += "# " +  FIR.makeURL(document.appform.elements, window.location.href + "\n");
        resultString += "out=fir(in";
        for (i=0; i<h.length; i++) {
            // Round the coefficients because BasicDSP has
            // a limit on line length.
            resultString += "," + FIR.roundNumber(h[i], 7);
        };
        resultString += ")\n";
        return resultString;
    },
    formatRaw: function(h) {
        var resultString = "";
        var i;
        for (i=0; i<h.length; i++) {
            resultString += h[i] + "\n";
        };
        return resultString;
    },
    // The Basic code that I started with for a band pass filter
    // is here:
    // http://www.dspguide.com/ch16/4.htm
    // (scroll down past the low-pass code
    // to see the band-pass code)
    // It 
    //   1. starts with low-pass filters for the low and high
    //       cutoff frequencies, 
    //   2. inverts the spectrum of the high cutoff frequency low-pass filter
    //      to make it a high-pass filter, 
    //   3. adds the low-pass and high kernels together to make
    //       a band-reject filter and
    //   4. inverts the spectrum of the band-reject filter to turn
    //       it into a band *pass* filter.
    // See also:
    //   Chapter 14: Introduction to Digital Filters
    //   http://www.dspguide.com/ch14/5.htm
    calculateFilter: function(samplerate, filterorder, 
            cutoff1, cutoff2, filtertype) {
        if (filtertype[0].checked) {      // low-pass
            FIR.filtertype = "low-pass";
            return FIR.calculateLowpass(samplerate, filterorder, cutoff1);
        }
        else if (filtertype[1].checked) { // high-pass
            FIR.filtertype = "high-pass";
            return FIR.calculateHighpass(samplerate, filterorder, cutoff1);
        }
        else if (filtertype[2].checked) { // band-reject
            FIR.filtertype = "band-reject";
            return FIR.calculateBandreject(samplerate, filterorder,
                    cutoff1, cutoff2);
        }
        else {                            // band-pass
            FIR.filtertype = "band-pass";
            return FIR.invertSpectrum(FIR.calculateBandreject(samplerate,
                    filterorder, cutoff1, cutoff2));
        };
    },
    graphdata: function(d, where) {
        var i;
        var g = [];
        for (i=0; i<d.length; i++) {
            g.push([i, d[i]]);
        };
        $(function () {
            $.plot($(where), [ g ]);
        });
    },
    graphfr: function(samplerate, d, where) {
        var i;
        var g = [];
        for (i=0; i<d.length; i++) {
            g.push([samplerate*(i/(d.length-1)/2),
                    // convert the frequency response to decibels
                    // This is explained at the bottom of
                    // http://www.dspguide.com/ch14/1.htm
                    20* Math.LOG10E*Math.log(d[i])]);
        };
        $(function () {
            $.plot($(where), [ { data: g, label: "dB" } ], { yaxis: { min: -100 } });
        });
    },
    makeURL: function(elements, baseURL) {
        var outURL = baseURL.split("?")[0];
        var i;
        for (i=0; i<elements.length; i++) {
            if (elements[i].name=="filtertype") {
                if (elements[i].checked) {
                    outURL += (i==0 ? "?" : "&") + elements[i].name + "=" + elements[i].value;                };
            }
            else {
                outURL += (i==0 ? "?" : "&") + elements[i].name + "=" + elements[i].value;
            };
        };
        return outURL;
    },
    // roundNumber is from
    // http://www.mediacollege.com/internet/javascript/number/round.html
    roundNumber: function(rnum, rlength) { // Arguments: number to round, number of decimal places
        return Math.round(rnum*Math.pow(10,rlength))/Math.pow(10,rlength);
    },
    isDefined: function(object, variable)
    // from http://jehiah.cz/archive/javascript-isdefined-function
    {
        return (typeof(eval(object)[variable]) == "undefined") ? false: true;
    },
    init: function(url, form) {
        var i;
        var j;
        var k;
        var v;
        var args;
        var argsx;
        if (url.indexOf("?") != -1) {
            args = url.split("?")[1];
            argsx = args.split("&");
            for (i=0; i<argsx.length; i++) {
                k = argsx[i].split("=")[0];
                v = argsx[i].split("=")[1];
                if ( FIR.isDefined(form,k) ) {
                    if (k=="filtertype") {
                        for (j=0; j<4; j++) {
                            if (form.filtertype[i+j-2].value==v) {
                                form.filtertype[i+j-2].checked=true;
                            }
                            else {
                                form.filtertype[i+j-2].checked=false;                            
                            };
                        };
                    }
                    else {
                        form[k].value = v;
                    }
                };
            };
        };
    },
    go: function() {
        var h = FIR.calculateFilter(document.appform.samplerate.value,
                          document.appform.filterorder.value,
                          document.appform.cutoff.value,
                          document.appform.highcutoff.value,
                          document.appform.filtertype);
        // We're doing two things in the next statement below:
        //   1. We convert the filter's kernel (impulse response) to its
        //       frequency response by applying the Discrete Fourier Transform (DFT)
        //   2. But *first* we add zeros to the the kernel to increase the resolution
        //       of the resulting frequency response. If I don't do this I don't see the
        //       ripple in the stop band.
        // Somewhere in The Scientist and Engineer's Guide to Digital Signal Processing
        // there is an explanation of this, but I can't find it at the moment.
        var fr = FIR.DFT(h.concat(FIR.arrayOfZeros(h.length*1)));
        document.outform.basicdspcoefficients.value = 
            FIR.formatBasicDSP(document.appform.samplerate.value,
                          document.appform.filterorder.value,
                          document.appform.cutoff.value,
                          document.appform.highcutoff.value,
                          document.appform.filtertype, h);
        document.outform.rawcoefficients.value =
            FIR.formatRaw(h);
        document.outform.frequencyresponse.value =
            FIR.formatRaw(fr.mag);
        document.outform.url.value =
            FIR.makeURL(document.appform.elements, window.location.href);
        FIR.graphdata(h,"#impulseresponse");
        FIR.graphfr(document.appform.samplerate.value, fr.mag, "#frequencyresponse");
    }
};
</script>
</head>
<body>
<h1>Filter Design - Calculate and display the impulse response (kernel) and frequency response for a windowed-sinc FIR filter - firkernel.htm</h1>
<form name="appform">
<label for="samplerate">Sample rate: </label> 
<input type="text" name="samplerate" size="10" onChange="FIR.go()" value="8000">
<br />
<label for="filterorder">Filter order: </label> 
<input type="text" name="filterorder" size="10" onChange="FIR.go()" value="100">
<br />
<label for="filtertype">Filter type: </label>
<input type="radio" name="filtertype" onChange="FIR.go()" value="lowpass">Low-pass
<input type="radio" name="filtertype" onChange="FIR.go()" value="highpass">High-pass
<input type="radio" name="filtertype" onChange="FIR.go()" value="bandreject">Band-reject
<input type="radio" name="filtertype" onChange="FIR.go()" value="bandpass" checked>Band-pass
<br />
<label for="cutofffrequency">Cutoff frequency: </label> 
<input type="text" name="cutoff" size="10" onChange="FIR.go()" value="1000"> (Hz)
<br />
<label for="highcutofffrequency">Cutoff frequency: </label> 
<input type="text" name="highcutoff" size="10" onChange="FIR.go()" value="2000"> (Hz - Band-pass & Band-reject)
<p>
</form>
<form name="outform">
<input type="button" onclick="FIR.go();" value="Calcluate coefficients">
</p>
<h3>Kernel (Impulse Response)</h3>
    <div id="impulseresponse" style="width:400px;height:200px;"></div>
<h3>Frequency Response</h3>
    <div id="frequencyresponse" style="width:400px;height:200px;"></div>

<br />
<label for="basicdspcoefficients">BasicDSP fir() Coefficients: </label>(
<a href="http://wwwhome.cs.utwente.nl/~ptdeboer/ham/basicdsp/">http://wwwhome.cs.utwente.nl/~ptdeboer/ham/basicdsp/</a> )
<br />
<textarea name="basicdspcoefficients" cols="80" rows="12"></textarea>
<br />
<label for="rawcoefficients">Kernel (Impulse Response) Raw Coefficients: </label>
<br />
<textarea name="rawcoefficients" cols="40" rows="12"></textarea>
<br />
<label for="freqencyresponse">Frequency response: </label>
<br />
<textarea name="frequencyresponse" cols="40" rows="12"></textarea>
<br />
<label for="url">url: </label>
<br />
<textarea name="url" cols="40" rows="6"></textarea>
<br />
<input type="button" onclick="window.location.href = document.outform.url.value;" value="Go">
</form>

<script id="source" language="javascript" type="text/javascript">
$(function () {
    var d1 = [];
    for (var i = 0; i < 14; i += 0.5)
        d1.push([i, Math.sin(i)]);

    var d2 = [[0, 3], [4, 8], [8, 5], [9, 13]];

    // a null signifies separate line segments
    var d3 = [[0, 12], [7, 12], null, [7, 2.5], [12, 2.5]];
    
    $.plot($("#impulseresponse"), [ d1 ]);
    $.plot($("#frequencyresponse"), [ d1 ]);
});
FIR.init(window.location.href, document.appform)
FIR.go();
</script>
<h2>Why?</h2>
Filters are basic to electronics, as bearings are to mechanical engineering. I created this page for several reasons:
<ol>
  <li>to better understand how filters work. Much of what I understand about filters I learned from <a id="dspbook" href="http://www.dspguide.com/">The Scientist and Engineer's Guide to Digital Signal Processing</a>. I recommend it.</li>
  <li>to create FIR (Finite Impulse Response) filter kernels that I can copy and paste 
into <a href="http://wwwhome.cs.utwente.nl/~ptdeboer/ham/basicdsp/">BasicDSP</a> (which has a FIR function), 
"A program for experimenting with simple audio DSP algorithms". 
I found <a href="http://home.tiscali.nl/curious_about/PE1OIT/basicdsp/index.html">BasicDSP</a> by <a href="http://www.google.com/search?rlz=1C1GGLS_enUS304US304&sourceid=chrome&ie=UTF-8&q=PA3FWM">Googling PA3FWM</a>, the author of <a href="http://www.websdr.org/">WebSDR</a>. PA3FWM (Pieter-Tjerk de Boer) is
also one of the two authors of BasicDSP.
(WebSDR is "a Software-Defined Radio receiver connected to the internet, 
allowing many listeners to listen and tune it simultaneously." WebSDR illuminates the use of
a filter to extract a signal from the ether.</li>
  <li>I have a long-standing interest (dare I say obsession) with 
<a href="http://en.wikipedia.org/wiki/Dual-tone_multi-frequency">DTMF</a> - Dual-Tone Multi-Frequency signally. I've used 
filters I designed with this page in BasicDSP to create a system that separates the low and high frequency bands used in DTMF.</li>
</ol>
<h2>Notes</h2>
<ul>
<li>When you design a filter with this page you get (in addition to the kernel and frequency response) a URL that will get your design back.</li>
<li>The JavaScript in this page is heavily commented, with references to <a href="#dspbook">the DSP book</a>, to explain how it works. The
source is available as a <a href="http://docs.google.com/Doc?id=accnfh3pfb3_361g77gvbdn&hl=en">Google Docs document</a>. (The browser on my Palm doesn't have View/Source.) - <a href="http://docs.google.com/Doc?id=accnfh3pfb3_361g77gvbdn&hl=en">http://docs.google.com/Doc?id=accnfh3pfb3_361g77gvbdn&hl=en</a></li>
<li>Other pages I've found helpful in understanding filters:
  <ul>
  <li><a href="http://www.falstad.com/">Paul Falstad's</a> <a href="http://www.falstad.com/dfilter/">Digital Filters Applet</a> - 
It allows you to play with many different types of filters. It graphically demonstrates the relationship between the impulse response
and frequency response, but it doesn't make available the impulse response (kernel) coefficients. - <a href="http://www.falstad.com/dfilter/">http://www.falstad.com/dfilter/</a></li>
  <li>MAXIM's <a href="http://www.maxim-ic.com/appnotes.cfm/appnote_number/733">A Filter Primer</a> - More than any other source it helped me understand poles and zeros and their effect on ripple in Buttterworth
and Chebychev filter design. <a href="http://www.maxim-ic.com/appnotes.cfm/appnote_number/733">http://www.maxim-ic.com/appnotes.cfm/appnote_number/733</a></li>
  <li><a href="http://quickfiltertech.com/">Quick Filter Technologies</a> offers a series of FIR intergrated circuits. Since their entire business is based on FIR filters their
website has a lot of information about them, including a 
<a href="http://www.quickfiltertech.com/html/content_page.php?content_id=33&">free Microsoft 
Windows FIR filter design download</a> that I discovered after
I started working on this page. <a href="http://quickfiltertech.com/">http://quickfiltertech.com/</a> | 
<a href="http://www.quickfiltertech.com/html/content_page.php?content_id=33&">http://www.quickfiltertech.com/html/content_page.php?content_id=33&</a></li>
<li><a href="http://del.icio.us/mshook/filter">http://del.icio.us/mshook/filter</a>
or <a href="http://mshook.googlepages.com/d4m.htm?/mshook/filter">http://mshook.googlepages.com/d4m.htm?/mshook/filter</a></li>
  </ul>
<li>I used <a href="http://code.google.com/p/flot/">FLOT</a> to create the graphs. I found it made creating attractive graphs very easy. FLOT, in turn, uses <a href="http://jquery.com/">jQuery</a>.</li>
<li>Comments or questions: <a href="mailto:mshook@gmail.com">mshook@gmail.com</a></li>
</ul>
<p>Tags: dsp filter javascript fir</p>
</body>
</html>




